Paradoxes of Subdiffusive Infiltration in Disordered Systems 
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Infiltration of diffusing particles from one material to another where the diffusion mechanism is 
either normal or anomalous is a widely observed phenomena. When the diffusion is anomalous we 
find interesting behaviors: diffusion may lead to an averaged net drift (x) from one material to 
another even if all particles eventually flow in the opposite direction, or may lead to a flow without 
drift. Starting with an underlying continuous time random walk model we solve diffusion equations 
describing this problem. Similar drift against flow is found in the quenched trap model. We argue 
that such a behavior is a general feature of diffusion in disordered systems. 
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Infiltration of diffusing particles from one material to 
another is a widely investigated process in many fields 
of physics. In recent years much focus was diverted to 
the problem when the diffusion in one or in both mate- 
rials is anomalous, namely (a;^) ^ with a 7^ 1 [J 01 ■ 
Among many examples where this behavior is important 
are infiltration of water into porous soil [3] , contaminant 
diffusion Q, moisture ingress in zeolites [5] or in fired 
clay ceramics f6| , diffusion of sugar through a membrane 
in a gel solvent [71, and polymer translocation through a 
membrane pore Infiltration is also important in bi- 
ologically motivated experiments. For example proteins 
diffusion is anomalous in the cell and normal in the exte- 
rior, compartments on membranes indicate that diffusion 
of proteins is taking place between different regions with 
varied diffusion mechanisms morphogens are subd- 
iffusing in extracellular environment where the diffusive 
properties changes abruptly in space [lo|. 

Consider unbiased diffusion in one-dimension where 
one type of diffusion takes place in a; < and another 
in a; > 0. The infiltration of particles from one material 
to another may lead to an averaged net drift {x). We 
show that for subdiffusion the fiow of particles may be 
in the opposite direction of the drift. Even more sur- 
prisingly, we find situations when asymptotically all the 
particles are in one sample but the average drift (x) is 
oppositely directed. This is a paradoxical behavior in 
the following sense: let t) be the normalized prob- 
ability density function (PDF) of finding a particle at 
time t in (x, x -\- dx), P{x, t)dx = 1. One can argue 
rather generally that if lim(_j.oo /o°° P{x^ t)dx = 1, i.e. all 
particles are in x > 0, then \m\t-^oo ^'^^P{x,t)dx = 
which implies that P(x, t) = for x < (since P{x, t) > 
0) and hence the drift should be positively directed 
limf_>.oo (x) — xP{x, t)dx = xP{x, t)dx > 0. 

Paradoxically in some cases of subdiffusion we find the 
opposite behavior limt_>.oo (x) < (even though all par- 
ticles are eventually in x > 0). Similarly, in some cases 
limt_>.oo (x) = 0, but in the long time limit all parti- 
cles accumulate in one sample, for example in x > 0, 
Ymit^oo P{x,t)dx = 1. The solution of these para- 
doxes is given in this Letter, as well as a derivation of (x) 



and Jq°° Pix, t)dx. 

Model 1: Fractional diffusion equations 0, — 
Consider the semi-infinite region x < which has subd- 
iffusive dynamics with exponent < a~ < 1 and diffu- 
sion constant K~ whose units are mt^/sec" . Similarly 
for the domain x > the exponent < < 1 and 



K+ 



mt^/sec" govern the dynamics. Subdiffusive pro- 



cesses are described by fractional diffusion equations [1^ 
(see also [3) 



dP{x,t) 

di 
dP{x,t) 

dt 



= oDl- 



dx"^ 
dx 



P(x,t), X < 0, 



X > 0, 



(1) 



where the Riemann-Liouville operator is defined as |13| 
oDl-°'P{x,t) = r-i(a) d/dt J*dt'P{x,t') [t - t'Y'\ 
The fractional diffusion equation Eq. ([Ij with — a+ — 
a and = = K yields for particles starting on 
the origin (x^) = 2Kt°'/r{l + a). For a" = a+ = 
1 this equation reduces to standard diffusion equation. 
Without the boundary conditions (soon to be derived), 
Eq. (m is nearly useless. The underlying random walk 
model we consider is the continuous time random walk 
(CTRW) [H, i. El] which is now specified. 

Model 2: CTRW. — Consider a jump process on a dis- 
crete lattice with the lattice spacing a. For lattice points 
X < a particle has the probability 1/2 to jump to one of 
its nearest neighbors. Waiting times on each lattice point 
are independent identically distributed random variables 
with a common PDF ip~ (r) . For x > a similar unbiased 
random walk takes place with a waiting time PDF ip^^r). 
On the lattice point x = (the boundary) a particle has 
the probability to jump right g+ or left q~ = 1 — q'^ [TBj 
and the waiting times are exponentially distributed with 
a rate Rq. Such biased interface is due for example to 
a difference of chemical potentials between the two sam- 
ples [Tcj . Thus, a particle starting on the origin will jump 
say to the right (with prob. q~^) after waiting an average 
time l/i?Oi then on the lattice point x = a, it will wait 
for time r drawn from i/)"'"(r), and then with probability 
1/2 will jump to the left or right. For subdiffusion the 
waiting times have power law PDFs ^/'"(t) cx t~(^+" ^ 
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FIG. 1: (color online). The drift, (a;), and occupation frac- 
tion in X < 0, 'P~{t), calculated numerically for the CTRW 
model (open circles) with a"*" = 0.75, = 0.138, a~ = 0.3, 
K~ — 0.385, g"*" = 0.7 and for the quenched trap model (filled 
circles) with 7+ — 0.9, 7" — 0.3 and q"*" — 0.7 (averaged over 
20 realizations of disorder). Dashed lines represent long time 
asymptotic Eqs. ([5]) and (|16|) . The dashed-dotted line given 
by Eq. psp describes how V~ approach its limit —J- 1. 
Notice that all particles flow to the left (P" —J- 1), however 
{x{t)) > namely particles drift to the right. 



and iP^{t) (X t^('^+°'^\ a.s t 00. More specifically, 
using standard Tauberian theorem the Laplace trans- 
form T — > s of the waiting time PDFs behave like 
ip-{s) ~ 1 - 3-3°"' , il)+{s) ~ 1 - B+s°'^ when s ^ 
corresponding to r — > 00 AH along this work we de- 
note the Laplace transform by the variable in the paren- 
theses /(s) = dt e"''*/(t). The generalized diffusion 
constants are given by — limy2_j.o^s-_j.o a^/2i?^ and 
= limjj2_>o.s+^o fl^/2-B+ [17|- Our results are not 
changed if on a; = the waiting times are power law 
distributed like ip~ or instead of exponential. Large 
number of applications of the CTRW model are discussed 
in BiH. 

The drift (x). — Using the CTRW approach we now 
calculate the drift. The position of a particle is a; = 
X]t=o 'where 6xi is the ith displacement and N is the 
random number of steps. Since the motion is unbiased 
in domains a; < and a; > 0, we have 



{x{t)) ^a{q+ -q-){nS))^ 



(2) 



where {nz{t)) is the average number of times a parti- 
cle visited the origin. We define a three state process 
^(t) = if the particle is on the origin, £^{t) = if 
the particle is in x > and £^{t) = —1 if the particle 
is in a: < 0. In the long time limit the number of vis- 
its to the origin is independent of i?o since the average 
waiting times in state -f and — are infinite. The wait- 
ing times in states -|- and — are the first passage times 



|18| from a; = a to a; = and from —a to 0, respec- 
tively. These first passage times in the continuum limit 
are one sided Levy distributions whose long time (small 
s) Laplace transforms are [l^ 4>~{s) ^ 1 — as" /^/V 
for a; < and similarly </>"'' (s) ~ 1 — as"^ /'^ / K+ for 
X > 0. The Laplace transform of the probability to have 
exactly transitions to state $^{t) = is easily found 
using the Laplace transform convolution theorem (20l | 



s 

where ^(s) = q^(j)^{s) + (7+0+ (s). From Eq. (jS]) 

0(s) 



s(l-^(s)) 



(3) 



(4) 



Using the small s expansion of (rtz) and Eq. ^ we obtain 
for < 



r(l + a-/2) 



(5) 



which agrees well with simulation in Fig. [TJ Similar ex- 
pression is found for a"*" < a~ . The sign of the drift, i.e. 
its directionality, is determined by the sign of 17+ — q~ , 
and (a:) = if g"*" = q~ ■ Eq. ([5]) shows that in the long 
time limit the drift depends only on one diffusion con- 
stant in sample (— ) and grows in time with the exponent 
of the slower medium. This is a surprising result: {x{t)) 
can be very far from the interface, deep in the faster 
sample x > 0, but still is independent of the properties 
of that region a+, K^. 

Boundary Conditions and Solution of Model 1). — 
Using the initial condition given by P(x, 0) — 5{x) the 
solution of Eq. ([!]) in Laplace space is given by 



ga+ /2 Igxp 



P(x,s) = C+(s)- 



/2 1 



+C-{s)- 




0{x) 



[i-e{x)], (6) 



where 9{x) is the step function. To find C+(s) and (s) 
we need two boundary conditions. The conservation of 
probability / dxP{x, s) — 1/s gives 



C-{s) + C+{s) 



(7) 



From Eq. ([7]) we get the first boundary condition which 
is simply the conservation of the probability current at 
the boundary 



J+{x = 0+,t) - J-{x = 0-,t) = ld{t) 



(8) 



D^-°' dP{x,t)/dx for X < 



where J-{x,t) = -R- nD]-"' 

and similarly for x > |2l| . To derive the second bound 
ary condition we calculate the first moment, (x(s)) = 
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FIG. 2: (color online). PDF of particle's position P{x,t) 
calculated numerically by CTRW model with q~ — 0.5, 
K- = 0.282, Q+ = 0.75, K+ = 0.138 and q+ = 0.7 at time 
t = 10^ . Dashed lines represent solution of fractional diffusion 
equations Eq. ((T| calculated by numerical Laplace inversion 
of Eqs. (|6I10|) . Notice the jump at the boundary. 



J dx X P{x, s), using Eq. ^ 



(xis)) 



1 

2^ 



(9) 

We require Eq. (|9]) to be equal to {x{s)) Eq. ([5]) calculated 
from the CTRW model. For a" < a+, Eqs. ([STfig)) yields 
when s — >■ 



2q^ 



K- 



C-{s) =2-C+(s). (10) 



Inverting Eqs. (|6ll0p to the time domain, analytical so- 
lutions are in excellent agreement with numerical simu- 
lations of the underlying CTRW model (see Fig. [5]). For 
the special case a~ = a^, we get 



C+is) = 



C-{s)=2-C+{s). (11) 



Using Eqs. (I6I10|) or pT|) . we derive the second boundary 
condition 

q+K-s-"'P{x = 0", s) q-R+s-"^ P{x ^ 0+, s), 

(12) 

which shows that generally the PDF at the boundary is 
not continuous, similar to the normal diffusion case [15l |. 
Such behavior is demonstrated in Fig. [21 Note that the 
boundary condition Eq. (jl2p has a form of convolution 
when ^ a~ . 

Occupation Fractions. — From Eqs. (j6ll0p occupation 
fractions, that is probabilities to be in x < or in 
a; > are V' {s) = J°^dx P{x,s) = C"(s)/(2s) and 



r+{s) = J^dx P{x,s) = C+(s)/(2s). In the long time 
limit equivalent to s — > 



1 



s(l + 7^-l(s)) 



, V-is) 



1 



(13) 



with 7^(s) = {q+VK^)/{q-\^) s("^-" Let us 

assume a~ < Since V^{s) ^ 1/s when s — >■ 0, then 
the probability to be in a; < is 



v-it)^i, t 



(14) 



which indicates that in the long time limit all particles 
flow to the region a; < 0, where the diffusion is slower (see 
Fig. [T]). Similar result can be obtained from the CTRW 
model [16.] 



r-{t) 




(15) 



Eq. (lisp is in excellent agreement with numerical simu- 
lations (see Fig. [T]). 

Paradox and its explanation. — (i) We observe an av- 
eraged net drift from one material to another even if 
all particles eventually flow in the opposite direction. If 
a~ < a"*" and q^ < q^ the drift {x) is positive (see Eq. 
([5]) and Fig. [1]). However, the particles are accumulating 
in X < 0, that \s V ^ I (see Eq. (HH), Fig. [1]). So, we 
have the drift directed opposite to the flux of the particles 
even if all particles are eventually in the slower sample. 
(ii) If a~ < but the boundary is unbiased q~ = g"*", 
in the long time limit the particles will be concentrated 
in the region x < Q. However, since q~ = q~^ there is no 
drift, (x) — 0. An explanation of these paradoxes is as 
follows: Although the region with smaller a will accumu- 
late more and more particles in the long time limit there 
will be always some particles in the opposite region where 
a is larger. These particles are moving more freely and 
travel far away from the interface which will compensate 
the accumulation of particles in the region with smaller a. 
In other words while = 1 — V~ — P[x, t)dx — >■ 0, 
X P{x,t)dx does not approach zero. While dynam- 
ics in faster domain {x > 0) is clearly important (since 
(x) > 0), note that (x) is independent of the diffusion 
properties of domain x > 0. Surprisingly, (x) Eq. ([5]) 
does not depend on a'^ and K'^ as mentioned. 

Model 3: Quenched trap model. — We proceed to show 
that effects discussed for the CTRW model and frac- 
tional diffusion equation are found also for systems with 
quenched disorder. Consider the quenched trap model 
where a particle is undergoing a one-dimensional ran- 
dom walk on a quenched random energy landscape on 
a lattice [H, [l^. On each lattice point a random en- 
ergy Ex is assigned, which is minus the energy of the 
particle on site x, so i^^r > is the depth of a trap on 
site X. The energies of the traps are independent identi- 
cally distributed random variables with a common PDF 
p{E) = {1/Tg)exp{—E/Tg). Once the energy at some 
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site X is defined it stays constant in time (quenched dis- 
order), which makes a difference with the corresponding 
"annealed" CTRW problem. The lattice is coupled to a 
heat bath with temperature T which leads to particles 
escape from site x and jumps to one of its nearest neigh- 
bors. The average time it takes the particle to escape 
from site x is given by Arrhenius law — exp^E^/T). 
A small change in leads to exponential change in t^. 
The PDF of the waiting times can be easily calculated 
i/,(r) = 7r-(i+'') with 7 = T/Tg and t > 1. For 7 < 1 
(T < Tg), all the moments of ^(t) diverge which leads to 
anomalous diffusion [l|, aging (23j , and nontrivial occu- 
pation times [13] ■ After averaging over different realiza- 
tions of disorder the evolution of the PDF and the mean 
squared displacement (x^) of the model is described by 
the subdiffusive exponent given by 27/(1 + 7) [1]. 

We define our composite quenched trap model such 
that we have 7^ in a; < and 7^^ in a; > and an interface 
located a,t x — 0. On the boundary a particle has the 
probability to jump right g"*" or left = 1 — and the 
waiting times are exponentially distributed. Numerical 
simulations reveal the behavior that was found for CTRW 
model and fractional diffusion equations Eq. Namely, 
for 7~ < 7+ we find the drift (C is a constant) 

(x) ^C{q+ -q-)t^'^^^+^'\ (16) 

which is positive for (7+ > g^, while almost all particles 
are found in a; < 0, V~ — !■ 1 as shown in Fig. [1] 

Remark. — Boundary conditions of fractional diffusion 
equations and CTRW models are non-trivial and have 



attracted previous interest [lU. We note that the 
solution of fractional equation Eqs. (|6ll0p must be 
used with care. While this solution gives the cor- 
rect asymptotic behavior of the occupation fraction 
T'^ — > 1 (when a" < a+) (see Eq. (|14p '). the correc- 
tion term within the fractional framework is (t) ^ 

1 - q+Vl<^/{q-VK+) (1 + "'-"^ ^ t""^^. This 

correction term is different from the exact CTRW result 
Eq. ([15]) (compare the Gamma functions). Thus, frac- 
tional equation works in the long time limit and already 
leading corrections to asymptotic solution show devia- 
tions from exact result. 

To summarize, we investigated infiltration in subdif- 
fusive systems. Using the CTRW model we derived the 
boundary conditions of the problem which allow analyt- 
ical solution of the fractional diffusion equations. Parti- 
cles flow to the slower medium while the direction of the 
averaged drift is determined by breaking of symmetry, 
Ql ^ qR in our model. This leads to interesting phenom- 
ena unique to anomalous diffusion: (i) all particles are 
found in one sample {V^ — >■ 1), but the drift is oppo- 
sitely directed {{x) > 0), (ii) drift does not depend on 
properties of fast medium {{x) is independent of a"*", 
even though (a;) might be located deep in that medium). 
We observe similar behavior for the composite quenched 
trap model which points out to a broader generality of 
our results. 
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